from scipy.stats import norm
import numpy as np
import math
import matplotlib.pyplot as plt
def deg2rad(deg):
    return deg * math.pi / 180

# 计算目标当前的横向位置分布
def get_x_distribution(D, delta_D, phi, delta_phi, V, delta_V, Vt, delta_Vt, X, delta_beta, delta_C):
    # 根据当前的参数计算 phi 值
    k = (V - delta_V) / (Vt - delta_Vt)
    phi = math.asin(k * math.sin(X))
    # 计算 gamma 值，用于后面计算分母
    sign = 1 if beta-C >= 0 else -1
    X += delta_beta - delta_C
    gamma = X + phi - delta_phi + delta_C + sign * math.pi
    # 计算分母并避免分母为零或接近零的情况
    denominator = np.sin(X + phi - delta_phi + delta_C)
    if abs(denominator) < 1e-6:
        denominator = 1e-6
    # 计算 mu 和 sigma
    mu = (D - delta_D) * (np.sin(phi - delta_phi + delta_beta) - k * np.sin(gamma)) / denominator
    sigma = 5 * D
    # 返回横向位置的概率密度函数
    return norm(loc=mu, scale=sigma)

# 计算击中概率
def get_hit_prob(D, delta_D, phi, delta_phi, V, delta_V, Vt, delta_Vt, l, beta, delta_beta, C, delta_C):
    hit_prob = 0
    # 设置舷角角度的步长
    step = 0.1
    angles = []
    probabilities = []
    for i in range(0, 181, int(1/step)):
        # 将角度转化为弧度
        X = deg2rad(i)
        # 计算当前舷角下的击中概率，并累加到结果中
        x_dist = get_x_distribution(D, delta_D, phi, delta_phi, V, delta_V, Vt, delta_Vt, X, delta_beta, delta_C)
        prob = x_dist.cdf(l/2) - x_dist.cdf(-l/2)
        hit_prob += math.cos(X) * prob * step
        # 输出当前舷角角度和对应的击中概率
        # print("舷角角度: %.2f度，击中概率: %.4f" % (i, hit_prob))
        angles.append(i)
        probabilities.append(hit_prob)
    return angles, probabilities

# 示例参数
Ds = [5, 10, 15, 20]
delta_Ds = [np.random.normal(scale=0.1*D) for D in Ds]
phi = deg2rad(10)
delta_phi = np.random.normal(scale=0.1)
V = 15
delta_V = np.random.normal(scale=0.1)
Vt = 50
delta_Vt = np.random.normal(scale=0.1)
l = 150
beta = deg2rad(60)
delta_beta = np.random.normal(scale=0.1)
C = deg2rad(30)
delta_C = np.random.normal(scale=0.1)

# 计算并输出每个舷角的击中概率
plt.figure(figsize=(8, 6))
colors = ['b', 'g', 'r', 'c']
lines = ['-', '--', ':', '-.']
for i, (D, delta_D) in enumerate(zip(Ds, delta_Ds)):
    angles, probabilities = get_hit_prob(D, delta_D, phi, delta_phi, V, delta_V, Vt, delta_Vt, l, beta, delta_beta, C, delta_C)
    # 对概率小于等于 0 的点进行处理
    probabilities = [prob if prob > 0 else 0 for prob in probabilities]
    # 绘制当前 D 值下的曲线
    plt.plot(angles, probabilities, color=colors[i], linestyle=lines[i], label=f"D={D}")
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体为中文黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
plt.xlim([0, 180])
plt.ylim([0, 1])
plt.xlabel("舷角角度 (度)")
plt.ylabel("击中概率")
plt.title("击中概率与舷角的关系图")
plt.legend()
plt.show()
